!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   FEMilaro -- Finite Element Method toolkit
!!
!! FEMilaro is free software; you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 3 of the License, or
!! (at your option) any later version.
!!
!! FEMilaro is distributed in the hope that it will be useful, but
!! WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
!! General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with FEMilaro; If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>

!> Generic grid, \f$\underline{x},\underline{v}\in\mathbb{R}^d\f$
!!
!! \n
!!
!! In this module we consider a tensor product, phase space, Cartesian
!! grid, composed by elements and sides. In general we have to deal
!! with the following objects:
!! <ul>
!!  <li> elements in the space or velocity grid:
!!  \f$K_x\in\mathcal{T}_x\f$ and \f$K_v\in\mathcal{T}_v\f$ 
!!  <li> elements in the phase space grid: \f$\mathcal{T}\ni
!!  K=K_x\times K_v\f$
!!  <li> sides in the space or velocity grid:
!!  \f$s_x=K_x\cap K^\prime_x \in\mathcal{E}_x\f$ and
!!  \f$s_v=K_v\cap K^\prime_v \in\mathcal{E}_v\f$
!!  <li> sides in the phase space grid: 
!!  \f$s_{v,x}=s_x\times K_v\f$ and
!!  \f$s_{x,v}=K_x\times s_v\f$.
!! </ul>
!! These entities are related related by the following bottom-up set
!! operations:
!! \f{displaymath}{
!!  \begin{array}{rcl}
!!   \cap & : &
!!     \left\{\begin{array}{ccl}
!!      K_x\cap K^\prime_x &\mapsto & s_x \\
!!      K_v\cap K^\prime_v &\mapsto & s_v 
!!     \end{array}\right. \\[4mm]
!!   \times & : &
!!     \left\{\begin{array}{ccl}
!!      K_x\times K_v &\mapsto & K \\
!!      s_x\times K_v &\mapsto & s_{v,x} \\
!!      K_x\times s_v &\mapsto & s_{x,v} \\
!!     \end{array}\right. \\
!!  \end{array}
!! \f}
!! as well as by the following top-bottom ones:
!! \f{displaymath}{
!!  \begin{array}{rcl}
!!   \partial_x & : &
!!    \begin{array}{ccl}
!!     \partial_x K &\mapsto&
!!      \partial K_x\times K_v = \bigcup_{i=1}^{2d} s_{v,x}^i
!!    \end{array} \\[3mm]
!!   \partial_v & : &
!!    \begin{array}{ccl}
!!     \partial_v K &\mapsto&
!!      K_x\times \partial K_v = \bigcup_{i=1}^{2d} s_{x,v}^i
!!    \end{array} \\[3mm]
!!   \partial & : &
!!    \left\{\begin{array}{ccl}
!!     \partial K &\mapsto& \begin{array}{c}
!!        \partial_x K\cup\partial_v K \\
!!        \partial K_x\times K_v\cup K_x\times \partial K_v \\
!!        \bigcup_{i=1}^{2d} s_{v,x}^i
!!      \cup \bigcup_{i=1}^{2d} s_{x,v}^i
!!      \end{array} \\[6mm]
!!     \partial K_x &\mapsto& \bigcup_{i=1}^{2d} s_x^i \\[2mm]
!!     \partial K_v &\mapsto& \bigcup_{i=1}^{2d} s_v^i.
!!    \end{array}\right.
!!  \end{array}
!! \f}
!!
!! \section Periodicity
!!
!! Concerning periodic boundary conditions, sides are not repeated, so
!! there are as many sides in each direction as elements, the last
!! element being connected to the first side along each direction.
module mod_tps_phs_grid

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_fu_manager, only: &
   mod_fu_manager_initialized, &
   new_file_unit

 use mod_octave_io, only: &
   mod_octave_io_initialized, &
   write_octave,   &
   read_octave,    &
   read_octave_al, &
   locate_var

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_tps_phs_grid_constructor, &
   mod_tps_phs_grid_destructor,  &
   mod_tps_phs_grid_initialized, &
   t_real_lists, &
   t_tps_phs_grid, new_tps_phs_grid, clear, &
   t_tps_grid, &
   locate_point, &
   write_octave, read_octave_al

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! public members

 ! pointer arrays
 type t_real_lists
  real(wp), allocatable :: x(:)
 end type t_real_lists

 type p_t_s1
  type(t_s1), pointer :: p => null()
 end type p_t_s1
 type p_t_e1
  type(t_e1), pointer :: p => null()
 end type p_t_e1
 type p_t_s
  type(t_s),  pointer :: p => null()
 end type p_t_s
 type p_t_e
  type(t_e),  pointer :: p => null()
 end type p_t_e

 !> one-field side
 !!
 !! Each side is defined by one index and one multiindex: the index
 !! denotes the normal direction (must be between 1 and \c d), the
 !! multiindex is used in the same way as for elements. In practice,
 !! (n,i) identify the first side of element i along the direction n.
 type t_s1
  integer :: d !< dimension
  !> normal direction
  integer :: n
  !> side index
  integer, allocatable :: i(:)
  !> box coordinates (2,d): start,end in each direction (the two
  !! values coincide in the \c n direction)
  real(wp), allocatable :: box(:,:)
  !> box width (d), 1 in the \c n direction
  real(wp), allocatable :: bw(:)
  !> side area
  real(wp) :: a
  !> connected elements
  integer :: ie(2)
  type(p_t_e1) :: e(2)
  integer :: o !< global counter
 end type t_s1

 !> one-field element
 !!
 !! An element is identified by a multiindex, defining its position in
 !! the corresponding one-field grid.
 type t_e1
  integer :: d !< dimension
  !> element multiindex (d)
  integer, allocatable :: i(:)
  !> box coordinates (2,d), start-end in each direction
  real(wp), allocatable :: box(:,:)
  !> box width
  real(wp), allocatable :: bw(:)
  !> element volume
  real(wp) :: vol
  !> connected sides (2,d), start-end in each direction
  integer, allocatable :: is(:,:)
  type(p_t_s1), allocatable :: s(:,:)
  !> sign indicator
  !!
  !! To each coordinate (evaluated at the barycenter) a sign indicator
  !! is associated as
  !! \f{displaymath}{
  !!  \iota_k = \left\{\begin{array}{ll}
  !!   1, & {\rm if}\quad x_{b_k}<0 \\
  !!   2, & {\rm if}\quad x_{b_k}\ge0
  !!  \end{array}\right.
  !! \f}
  !! and the element sign indicator is defined as
  !! \f{displaymath}{
  !!  \iota = 1 + \sum_{k=1}^d 2^{k-1}(\iota_k-1).
  !! \f}
  !! Notice that \f$\iota=1,\ldots\,2^d\f$.
  integer, allocatable :: si(:)
  integer :: o !< global counter
 end type t_e1

 !> phase space side
 !!
 !! This is a Cartesian product of a one-field side and a one-field
 !! element. In principle, all the properties can be recovered from
 !! such two entities, however it is convenient precomputing some of
 !! them.
 type t_s
  integer :: d !< dimension
  type(t_s1), pointer :: s1 => null() !< parent side
  type(t_e1), pointer :: e1 => null() !< parent element
  !> connected elements
  integer :: ie(2)
  type(p_t_e) :: e(2)
  integer :: o !< global counter
 end type t_s

 !> phase space element
 type t_e
  integer :: d !< dimension
  type(p_t_e1) :: e1(2) !< parent elements
  !> connected sides (start,end in each direction)
  type(p_t_s), allocatable :: sv_x(:,:)
  type(p_t_s), allocatable :: sx_v(:,:)
  integer :: o !< global counter
 end type t_e

 !> One-field grid (i.e. either space or velocity)
 type t_tps_grid
  integer :: d
  !> total number of elements
  integer :: ne
  !> number of elements in each direction
  integer, allocatable :: nen(:)
  !> total number of sides
  integer :: ns
  !> number of sides in each direction
  integer, allocatable :: nsn(:)
  !> sides
  type(t_s1), allocatable :: s(:)
  !> elements
  type(t_e1), allocatable :: e(:)
  !> private field to speed-up element indexing (see \c i2ie1)
  integer, private, allocatable :: esrc(:)
  integer, private, allocatable :: ssrc(:)
  integer, private :: ssrcn
 contains
  !> side index, given a pair (n,i)
  procedure, pass(grid) :: ni2is => ni2is1
  !> element index, given a multiindex
  procedure, pass(grid) :: i2ie  => i2ie1
 end type t_tps_grid

 !> Phase-space grid
 !!
 !! \note Since the complexity of the generic dimension is dealt with
 !! in \c t_tps_grid, this grid is much simpler and conceptually
 !! analogous to a two.dimensional grid, where the two dimensions are
 !! space and velocty. This fact is exploited to define the element
 !! storage order.
 type t_tps_phs_grid
  integer :: d
  type(t_tps_grid) :: gx !< parent grid (space)
  type(t_tps_grid) :: gv !< parent grid (velocity)
  !> total number of elements
  integer :: ne
  !> total number of sides: <tt>nsvx + nsxv</tt>
  integer :: ns
  !> number of sides in \c sv_x
  integer :: nsvx
  !> number of sides in \c sx_v
  integer :: nsxv
  !> elements
  type(t_e), allocatable :: e(:)
  !> \f$\partial_x K\f$
  type(t_s), allocatable :: sv_x(:)
  !> \f$\partial_v K\f$
  type(t_s), allocatable :: sx_v(:)
 contains
  !> side index, given a triplet (n,isx,iev)
  procedure, pass(grid) :: i2isvx => i2isvx
  !> side index, given a triplet (iex,n,isv)
  procedure, pass(grid) :: i2isxv => i2isxv
  !> element index, given a pair (iex,iev) of multiindexes
  procedure, pass(grid) :: i2ie  => i2ie
 end type t_tps_phs_grid

 ! private members

! Module variables

 ! public members
 logical, protected ::               &
   mod_tps_phs_grid_initialized = .false.
 ! private members
 character(len=*), parameter :: &
   this_mod_name = 'mod_tps_phs_grid'

 interface write_octave
   module procedure write_grid,  write_e,  write_s, &
                    write_grid1, write_e1, write_s1
 end interface write_octave

 interface read_octave_al
   module procedure read_rlist_al
 end interface read_octave_al

 interface clear
   module procedure clear_grid
 end interface

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_tps_phs_grid_constructor()
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.)   .or. &
       (mod_kinds_initialized.eqv..false.)      .or. &
       (mod_fu_manager_initialized.eqv..false.) .or. &
       (mod_octave_io_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_tps_phs_grid_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   mod_tps_phs_grid_initialized = .true.
 end subroutine mod_tps_phs_grid_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_tps_phs_grid_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_tps_phs_grid_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_tps_phs_grid_initialized = .false.
 end subroutine mod_tps_phs_grid_destructor

!-----------------------------------------------------------------------

 !> Define a new phase-space grid
 !!
 !! The grid is defined according to the coordinate vectors \c xx, \c
 !! vv, and periodic boundary conditions are assumed.
 !!
 !! The ordering of the phase-space elements and sides is completely
 !! controlled by \c i2ie and \c i2isvx, \c i2isxv: this function
 !! works for any chosen ordering.
 subroutine new_tps_phs_grid(grid,xx,vv)
  type(t_real_lists), intent(in) :: xx(:), vv(:)
  type(t_tps_phs_grid), intent(out), target :: grid

  integer :: ie, is, iex, iev, isx, isv, i
  type(t_e1), pointer :: e11, e12
  character(len=*), parameter :: &
    this_sub_name = 'new_tps_phs_grid'

   if(size(xx).ne.size(vv)) call error(this_sub_name,this_mod_name, &
      'The number of space and velocity dimensions must be the same')

   grid%d = size(xx)

   ! build the two one-field grids
   call new_grid1(grid%gx,xx)
   call new_grid1(grid%gv,vv)
   
   ! now the phase-space grid
   grid%ne = grid%gx%ne*grid%gv%ne
   grid%ns = grid%ne*(grid%gx%d+grid%gv%d) ! see new_grid1
   grid%nsvx = grid%gx%ns*grid%gv%ne
   grid%nsxv = grid%gx%ne*grid%gv%ns
   allocate( grid%e(grid%ne),grid%sv_x(grid%nsvx),grid%sx_v(grid%nsxv) )

   ! Fill the sides
   sv_x_side_do: do iev=1,grid%gv%ne
     do isx=1,grid%gx%ns
       is = grid%i2isvx( grid%gx%s(isx)%n , grid%gx%s(isx)%i , &
                         grid%gv%e(iev)%i )

       grid%sv_x(is)%d = grid%d
       grid%sv_x(is)%s1 => grid%gx%s(isx)
       grid%sv_x(is)%e1 => grid%gv%e(iev)
       grid%sv_x(is)%o = is

       ! To find the two elements, notice that they have to intersect
       ! along the x side, and share the same parent v element.
       grid%sv_x(is)%ie(1) = grid%i2ie(                  &
          grid%sv_x(is)%s1%e(1)%p%i , grid%sv_x(is)%e1%i )
       grid%sv_x(is)%ie(2) = grid%i2ie(                  &
          grid%sv_x(is)%s1%e(2)%p%i , grid%sv_x(is)%e1%i )
       grid%sv_x(is)%e(1)%p => grid%e( grid%sv_x(is)%ie(1) )
       grid%sv_x(is)%e(2)%p => grid%e( grid%sv_x(is)%ie(2) )
     enddo
   enddo sv_x_side_do
   sx_v_side_do: do isv=1,grid%gv%ns
     do iex=1,grid%gx%ne
       is = grid%i2isxv( grid%gx%e(iex)%i ,                    &
                         grid%gv%s(isv)%n , grid%gv%s(isv)%i )

       grid%sx_v(is)%d = grid%d
       grid%sx_v(is)%s1 => grid%gv%s(isv)
       grid%sx_v(is)%e1 => grid%gx%e(iex)
       grid%sx_v(is)%o = is

       ! To find the two elements, notice that they have to intersect
       ! along the v side, and share the same parent x element.
       grid%sx_v(is)%ie(1) = grid%i2ie(                  &
          grid%sx_v(is)%e1%i , grid%sx_v(is)%s1%e(1)%p%i )
       grid%sx_v(is)%ie(2) = grid%i2ie(                  &
          grid%sx_v(is)%e1%i , grid%sx_v(is)%s1%e(2)%p%i )
       grid%sx_v(is)%e(1)%p => grid%e( grid%sx_v(is)%ie(1) )
       grid%sx_v(is)%e(2)%p => grid%e( grid%sx_v(is)%ie(2) )
     enddo
   enddo sx_v_side_do

   ! Fill the elements
   do iev=1,grid%gv%ne
     do iex=1,grid%gx%ne
       ie = grid%i2ie( grid%gx%e(iex)%i , grid%gv%e(iev)%i )
       
       grid%e(ie)%d = grid%d
       grid%e(ie)%o = ie
       grid%e(ie)%e1(1)%p => grid%gx%e(iex)
       grid%e(ie)%e1(2)%p => grid%gv%e(iev)

       ! total number of connected sides: the sides of both parent
       ! elements become sides of the phase-space element

       ! sv_x: obtained from the sides of e1(1)
       e11 => grid%e(ie)%e1(1)%p
       e12 => grid%e(ie)%e1(2)%p
       allocate( grid%e(ie)%sv_x( 2 , size(e11%s,2) ) )
       do is=1,size(grid%e(ie)%sv_x,2)
         do i=1,2
           grid%e(ie)%sv_x(i,is)%p => grid%sv_x(                      &
             grid%i2isvx( e11%s(i,is)%p%n , e11%s(i,is)%p%i , e12%i ) )
         enddo
       enddo

       ! sx_v: obtained from the sides of e1(2)
       e11 => grid%e(ie)%e1(1)%p
       e12 => grid%e(ie)%e1(2)%p
       allocate( grid%e(ie)%sx_v( 2 , size(e12%s,2) ) )
       do is=1,size(grid%e(ie)%sx_v,2)
         do i=1,2
           grid%e(ie)%sx_v(i,is)%p => grid%sx_v(                      &
             grid%i2isxv( e11%i , e12%s(i,is)%p%n , e12%s(i,is)%p%i ) )
         enddo
       enddo

     enddo
   enddo

 end subroutine new_tps_phs_grid

!-----------------------------------------------------------------------
 
 pure subroutine clear_grid(grid)
  type(t_tps_phs_grid), intent(out) :: grid
  character(len=*), parameter :: &
    this_sub_name = 'clear_grid'
   
   ! nothing to do
 end subroutine clear_grid
 
!-----------------------------------------------------------------------

 subroutine new_grid1(grid,xx)
  type(t_real_lists), intent(in) :: xx(:)
  type(t_tps_grid), intent(out), target :: grid

  integer :: id, ie, is, e_cnt
  integer, allocatable :: midx(:), midxs(:)
  character(len=*), parameter :: &
    this_sub_name = 'new_grid1'

   grid%d = size(xx)

   ! Define the number of elements and sides
   allocate(grid%nen(grid%d),grid%nsn(grid%d))
   do id=1,grid%d
     grid%nen(id) = size(xx(id)%x)-1
     grid%nsn(id) = size(xx(id)%x)-1 ! periodicity: one el. -> one side
   enddo
   grid%ne = product(grid%nen)
   ! To count the sides, observe that (thanks to periodicity) to each
   ! element correspond d sides, one for each direction
   grid%ns = grid%ne*grid%d
   allocate(grid%s(grid%ns))
   allocate(grid%e(grid%ne))

   ! Set the private fields
   allocate(grid%esrc(grid%d-1),grid%ssrc(grid%d-1))
   do id=1,grid%d-1
     grid%esrc(id) = product(grid%nen(1:id))
     grid%ssrc(id) = product(grid%nsn(1:id))
   enddo
   grid%ssrcn = product(grid%nsn)

   ! Fill elements and sides: the trick here is to respect the storage
   ! order defined by i2ie, and for this reason all the element and
   ! side indexes are generated as multiindexes and converted with
   ! i2ie, i2is. Concerning the sides, we exploit here the fact that
   ! each side is seen exactly once as the "left" side of an element,
   ! for a given direction.
   allocate(midx(grid%d),midxs(grid%d)); midx = 1
   e_cnt = 0 ! element counter, to exit the following loop
   do
     e_cnt = e_cnt + 1
    if(e_cnt.gt.grid%ne) exit
     ie = grid%i2ie(midx)

     ! fill the element
     grid%e(ie)%d = grid%d
     allocate(grid%e(ie)%i(grid%d))
     grid%e(ie)%i = midx
     grid%e(ie)%o = ie
     allocate( grid%e(ie)%box(2,grid%d) , grid%e(ie)%bw(grid%d) )
     allocate( grid%e(ie)%si(grid%d) )
     do id=1,grid%d
       grid%e(ie)%box(:,id) = xx(id)%x( midx(id):midx(id)+1 )
       grid%e(ie)%bw (  id) = grid%e(ie)%box(2,id)-grid%e(ie)%box(1,id) 
       if( (grid%e(ie)%box(1,id)+0.5_wp*grid%e(ie)%bw(id)) &
          .lt.0.0_wp) then
         grid%e(ie)%si(id) = 1
       else
         grid%e(ie)%si(id) = 2
       endif
     enddo
     grid%e(ie)%vol = product( grid%e(ie)%bw )
     allocate( grid%e(ie)%is(2,grid%d) , grid%e(ie)%s(2,grid%d) )
     do id=1,grid%d
       ! the first side has the same multiindex as the element
       is = grid%ni2is(id,midx)
       grid%e(ie)%is(1,id)   =         is
       grid%e(ie)%s (1,id)%p => grid%s(is)
       ! the second side is obtained incrementing the id entry in the
       ! multiindex; notice that this may be a periodic side
       midxs = midx; midxs(id) = midxs(id)+1
       if(midxs(id).gt.grid%nsn(id)) midxs(id) = 1
       is = grid%ni2is(id,midxs)
       grid%e(ie)%is(2,id)   =         is
       grid%e(ie)%s (2,id)%p => grid%s(is)
     enddo

     ! fill the "left" side
     do id=1,grid%d
       is = grid%ni2is(id,midx) ! same as previous loop
       
       grid%s(is)%d = grid%d
       grid%s(is)%n = id
       allocate( grid%s(is)%i(grid%d) )
       grid%s(is)%i = midx
       grid%s(is)%o = is
       allocate( grid%s(is)%box(2,grid%d) , grid%s(is)%bw(grid%d) )
       ! the box definition is taken from the attached element, except
       ! for the id direction
       grid%s(is)%box = grid%e(ie)%box
       grid%s(is)%box(:,id) = xx(id)%x( midx(id) ) ! repeated value
       grid%s(is)%bw = grid%e(ie)%bw
       grid%s(is)%bw (  id) = 1.0_wp
       grid%s(is)%a = product( grid%s(is)%bw )
       ! the second element is already available
       grid%s(is)%ie(2) = ie
       grid%s(is)%e (2)%p => grid%e(ie)
       ! the first element can be periodic
       midxs = midx; midxs(id) = midxs(id)-1
       if(midxs(id).eq.0) midxs(id) = grid%nen(id)
       grid%s(is)%ie(1) = grid%i2ie(midxs)
       grid%s(is)%e (1)%p => grid%e( grid%s(is)%ie(1) )
     enddo
   
     ! increment the multiindex
     call midx_p1(midx,grid%nen)
   enddo
   deallocate(midx,midxs)

 contains

  ! increment the multiindex
  pure subroutine midx_p1(ii,nn)
   integer, intent(in) :: nn(:) ! shape in each dimension
   integer, intent(inout) :: ii(:)

   integer :: id

    do id=1,grid%d
      ii(id) = ii(id) + 1
      if(id.eq.grid%d) exit ! no bounds on the last index
      if(ii(id).gt.nn(id)) then
        ii(id) = 1 ! next dimension
      else
        exit
      endif
    enddo
    
  end subroutine midx_p1

 end subroutine new_grid1

!-----------------------------------------------------------------------

 !> Define the sequential side ordering
 !!
 !! The sides are ordered in column major order (like the elements),
 !! with an additional, most external index given by the direction n.
 pure function ni2is1(grid,n,i) result(o)
  class(t_tps_grid), target, intent(in) :: grid
  integer, intent(in) :: i(:), n
  integer :: o

  logical, parameter :: debug = .true.

   if(debug) then
     if(any(i.le.0).or.any(i.gt.grid%nsn).or. &
        (n.le.0).or.(n.gt.grid%d)) then
       o = -huge(o)/10
       return
     endif
   endif

   o = i(1) + sum((i(2:)-1)*grid%ssrc) + (n-1)*grid%ssrcn
 end function ni2is1

!-----------------------------------------------------------------------

 !> Define the sequential element ordering.
 !!
 !! The elements are ordered starting from the first dimension onward,
 !! according to a "column major" ordering. This mimics a \d
 !! dimensional fortran array with shape given by <tt>grid\%nen</tt>.
 pure function i2ie1(grid,i) result(o)
  class(t_tps_grid), target, intent(in) :: grid
  integer, intent(in) :: i(:)
  integer :: o

  logical, parameter :: debug = .false.

   if(debug) then
     if(any(i.le.0).or.any(i.gt.grid%nen)) then
       o = -huge(o)/10
       return
     endif
   endif
   
   o = i(1) + sum((i(2:)-1)*grid%esrc)
 end function i2ie1

!-----------------------------------------------------------------------

 !> Define the sequential side ordering.
 pure function i2isvx(grid,n,isx,iev) result(o)
  class(t_tps_phs_grid), target, intent(in) :: grid
  integer, intent(in) :: n, isx(:), iev(:)
  integer :: o

   ! space is the fastest index
   !o = grid%gx%ni2is(n,isx) + (grid%gv%i2ie(iev)-1)*grid%gx%ns
   ! velocity is the fastest index
   o = (grid%gx%ni2is(n,isx)-1)*grid%gv%ne + grid%gv%i2ie(iev)

 end function i2isvx

 pure function i2isxv(grid,iex,n,isv) result(o)
  class(t_tps_phs_grid), target, intent(in) :: grid
  integer, intent(in) :: iex(:), n, isv(:)
  integer :: o

   ! space is the fastest index
   !o = grid%gx%i2ie(iex) + (grid%gv%ni2is(n,isv)-1)*grid%gx%ne
   ! velocity is the fastest index
   o = (grid%gx%i2ie(iex)-1)*grid%gv%ns + grid%gv%ni2is(n,isv)

 end function i2isxv

!-----------------------------------------------------------------------

 !> Define the sequential element ordering.
 pure function i2ie(grid,ix,iv) result(o)
  class(t_tps_phs_grid), target, intent(in) :: grid
  integer, intent(in) :: ix(:), iv(:)
  integer :: o

   ! space is the fastest index
   !o = grid%gx%i2ie(ix) + (grid%gv%i2ie(iv)-1)*grid%gx%ne
   ! velocity is the fastest index
   o = (grid%gx%i2ie(ix)-1)*grid%gv%ne + grid%gv%i2ie(iv)

 end function i2ie

!-----------------------------------------------------------------------
 
 !> Locate the element a point belongs to
 !!
 !! This subroutine is similar to the corresponding one in \c
 !! mod_grid, except for the definition of the barycentric
 !! coordinates. Here we compute, for each dimension,
 !! \f{displaymath}{
 !!  \xi^k = \frac{x^k-x^k_l}{x^k_r-x^k_l}
 !! \f}
 !! where \f$x^k_l\f$ and \f$x^k_r\f$ are coordinates of the element
 !! box. The point \f$\underline{x}\f$ belongs to the element if
 !! \f$\xi^k\in[0\,,1]\f$ for each \f$k\f$. The sign of \f$\xi^k\f$
 !! also gives an indication about where to search next.
 !!
 !! \todo This function could be optimized.
 pure subroutine locate_point(ie,x,grid,ie0,ip)
  integer, allocatable, intent(out) :: ie(:) !< elements
  real(wp),         intent(in) :: x(:,:) !< points as column vectors
  type(t_tps_grid), intent(in) :: grid   !< grid
  integer,              intent(in), optional :: ie0(:) !< initial guess
  integer, allocatable, intent(out), optional :: ip(:) !< loc. points
 
  real(wp), parameter :: toll = 100.0_wp*epsilon(0.0_wp)
  logical :: search_completed, need_restart, found_restart
  logical, allocatable :: checked(:)
  integer :: np, i, je, ml(2)
  integer, allocatable :: wip(:), wie(:)
  real(wp) :: xi(grid%d,2)
  character(len=*), parameter :: &
    this_sub_name = 'locate_point'
 
   np = size(x,2)
   allocate( wip(np) , wie(np) ) ! work arrays
   allocate( checked(grid%ne) )

   np = 0 ! counter
   point_do: do i=1,size(x,2)

     checked = .false.

     if(present(ie0)) then
       je = ie0(i)
     else
       je = 1
     endif
     search_completed = .false.
     do
      if(search_completed) exit

       ! We check here all the coordinates: this is not very
       ! efficient because one could locate one coordinate at a time,
       ! using the regularity of the grid. However, this
       ! implementation is more robust in case one wants to
       ! generalize the grid structure.

       ! xi coordinate
       xi(:,1) = ( x(:,i)-grid%e(je)%box(1,:) )/grid%e(je)%bw(:)
       xi(:,2) = 1.0_wp - xi(:,1)
       if(minval(xi).ge.-toll) then ! found the element
         np = np+1
         wip(np) = i
         wie(np) = je
         search_completed = .true.
       else
         checked(je) = .true.
         ml = minloc(xi) ! row and column of the minimum
         ! Currently there is no element-to-element information: we
         ! need to pass through the sides.
         je = grid%s( grid%e(je)%is(ml(2),ml(1)) )%ie(ml(2))
         if(je.lt.1) then         ! boundary condition
           need_restart = .true.
         elseif(checked(je)) then ! already checked element
           need_restart = .true.
         else                     ! je is a good candidate
           need_restart = .false.
         endif
         if(need_restart) then
           call get_next(je,found_restart)
           if(.not.found_restart) then ! the point can not be located
             search_completed = .true.
           endif
         endif
       endif
     enddo
     
   enddo point_do

   allocate( ie(np) )
   ie = wie(1:np)
   if(present(ip)) then
     allocate( ip(np) )
     ip = wip(1:np)
   endif

   deallocate( wip , wie , checked )
 contains

  pure subroutine get_next(ie,found)
   integer, intent(out) :: ie
   logical, intent(out) :: found

    found = .false.
    do ie=1,grid%ne
      if(.not.checked(ie)) then
        found = .true.
        exit
      endif
    enddo
  end subroutine get_next

 end subroutine locate_point
 
!-----------------------------------------------------------------------

 subroutine write_grid(grid,var_name,fu)
 ! octave output for the complete grid
  integer, intent(in) :: fu
  type(t_tps_phs_grid), intent(in) :: grid
  character(len=*), intent(in) :: var_name
 
  character(len=*), parameter :: &
    this_sub_name = 'write_grid'
 
   write(fu,'(a,a)')    '# name: ',var_name
   write(fu,'(a)')      '# type: struct'
   write(fu,'(a)')      '# length: 10' ! number of fields

   ! field 01 : d
   write(fu,'(a)')      '# name: d'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(grid%d,'<cell-element>',fu)

   ! field 02 : gx
   write(fu,'(a)')      '# name: gx'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(grid%gx,'<cell-element>',fu)

   ! field 03 : gv
   write(fu,'(a)')      '# name: gv'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(grid%gv,'<cell-element>',fu)

   ! field 04 : ne
   write(fu,'(a)')      '# name: ne'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(grid%ne,'<cell-element>',fu)

   ! field 05 : ns
   write(fu,'(a)')      '# name: ns'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(grid%ns,'<cell-element>',fu)

   ! field 06 : nsvx
   write(fu,'(a)')      '# name: nsvx'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(grid%nsvx,'<cell-element>',fu)

   ! field 07 : nsxv
   write(fu,'(a)')      '# name: nsxv'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(grid%nsxv,'<cell-element>',fu)

   ! field 08 : e
   write(fu,'(a)')      '# name: e'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(grid%e,'<cell-element>',fu)

   ! field 09 : sv_x
   write(fu,'(a)')      '# name: sv_x'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(grid%sv_x,'<cell-element>',fu)

   ! field 10 : sx_v
   write(fu,'(a)')      '# name: sx_v'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(grid%sx_v,'<cell-element>',fu)

 end subroutine write_grid

!-----------------------------------------------------------------------

 subroutine write_s(s,var_name,fu)
  integer, intent(in) :: fu
  type(t_s), intent(in) :: s(:)
  character(len=*), intent(in) :: var_name
 
  integer :: is
  character(len=*), parameter :: &
    this_sub_name = 'write_s'

   write(fu,'(a,a)')    '# name: ',var_name
   write(fu,'(a)')      '# type: struct'
   write(fu,'(a)')      '# length: 4' ! number of fields

   ! field 01 : d
   write(fu,'(a)')      '# name: d'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(s)
   do is=1,size(s)
     call write_octave(s(is)%d,'<cell-element>',fu)
   enddo

   ! field 02 : s1
   write(fu,'(a)')      '# name: s1'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(s)
   do is=1,size(s)
     call write_octave(s(is)%s1%o,'<cell-element>',fu)
   enddo

   ! field 03 : e1
   write(fu,'(a)')      '# name: e1'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(s)
   do is=1,size(s)
     call write_octave(s(is)%e1%o,'<cell-element>',fu)
   enddo

   ! field 04 : ie
   write(fu,'(a)')      '# name: ie'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(s)
   do is=1,size(s)
     call write_octave(s(is)%ie,'r','<cell-element>',fu)
   enddo

 end subroutine write_s

!-----------------------------------------------------------------------

 subroutine write_e(e,var_name,fu)
  integer, intent(in) :: fu
  type(t_e), intent(in) :: e(:)
  character(len=*), intent(in) :: var_name
 
  integer :: ie, i, j
  character(len=*), parameter :: &
    this_sub_name = 'write_e'

   write(fu,'(a,a)')    '# name: ',var_name
   write(fu,'(a)')      '# type: struct'
   write(fu,'(a)')      '# length: 4' ! number of fields

   ! field 01 : d
   write(fu,'(a)')      '# name: d'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(e)
   do ie=1,size(e)
     call write_octave(e(ie)%d,'<cell-element>',fu)
   enddo

   ! field 02 : e1
   write(fu,'(a)')      '# name: e1'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(e)
   do ie=1,size(e)
     call write_octave( (/(e(ie)%e1(i)%p%o, i=1,2)/) , &
                                'r','<cell-element>',fu)
   enddo

   ! field 03 : sv_x
   write(fu,'(a)')      '# name: isv_x'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(e)
   do ie=1,size(e)
     call write_octave( reshape(                                      &
       (/(( e(ie)%sv_x(i,j)%p%o , i=1,2), j=1,size(e(ie)%sv_x,2))/) , &
       (/2,size(e(ie)%sv_x,2)/) ) , '<cell-element>',fu)
   enddo

   ! field 04 : sx_v
   write(fu,'(a)')      '# name: isx_v'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(e)
   do ie=1,size(e)
     call write_octave( reshape(                                      &
       (/(( e(ie)%sx_v(i,j)%p%o , i=1,2), j=1,size(e(ie)%sx_v,2))/) , &
       (/2,size(e(ie)%sx_v,2)/) ) , '<cell-element>',fu)
   enddo

 end subroutine write_e
 
!-----------------------------------------------------------------------

 subroutine write_grid1(grid,var_name,fu)
 ! octave output for the complete grid
  integer, intent(in) :: fu
  type(t_tps_grid), intent(in) :: grid
  character(len=*), intent(in) :: var_name
 
  character(len=*), parameter :: &
    this_sub_name = 'write_grid1'
 
   write(fu,'(a,a)')    '# name: ',var_name
   write(fu,'(a)')      '# type: struct'
   write(fu,'(a)')      '# length: 7' ! number of fields
 
   ! field 01 : d
   write(fu,'(a)')      '# name: d'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(grid%d,'<cell-element>',fu)

   ! field 02 : ne
   write(fu,'(a)')      '# name: ne'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(grid%ne,'<cell-element>',fu)

   ! field 03 : nen
   write(fu,'(a)')      '# name: nen'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(grid%nen,'r','<cell-element>',fu)

   ! field 04 : ns
   write(fu,'(a)')      '# name: ns'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(grid%ns,'<cell-element>',fu)

   ! field 05 : nsn
   write(fu,'(a)')      '# name: nsn'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(grid%nsn,'r','<cell-element>',fu)

   ! field 06 : s
   write(fu,'(a)')      '# name: s'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(grid%s,'<cell-element>',fu)

   ! field 07 : e
   write(fu,'(a)')      '# name: e'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(grid%e,'<cell-element>',fu)

 end subroutine write_grid1

!-----------------------------------------------------------------------

 subroutine write_s1(s,var_name,fu)
  integer, intent(in) :: fu
  type(t_s1), intent(in) :: s(:)
  character(len=*), intent(in) :: var_name
 
  integer :: is
  character(len=*), parameter :: &
    this_sub_name = 'write_s1'

   write(fu,'(a,a)')    '# name: ',var_name
   write(fu,'(a)')      '# type: struct'
   write(fu,'(a)')      '# length: 7' ! number of fields

   ! field 01 : d
   write(fu,'(a)')      '# name: d'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(s)
   do is=1,size(s)
     call write_octave(s(is)%d,'<cell-element>',fu)
   enddo

   ! field 02 : n
   write(fu,'(a)')      '# name: n'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(s)
   do is=1,size(s)
     call write_octave(s(is)%n,'<cell-element>',fu)
   enddo

   ! field 03 : i
   write(fu,'(a)')      '# name: i'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(s)
   do is=1,size(s)
     call write_octave(s(is)%i,'r','<cell-element>',fu)
   enddo

   ! field 04 : box
   write(fu,'(a)')      '# name: box'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(s)
   do is=1,size(s)
     call write_octave(s(is)%box,'<cell-element>',fu)
   enddo

   ! field 05 : bw
   write(fu,'(a)')      '# name: bw'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(s)
   do is=1,size(s)
     call write_octave(s(is)%bw,'r','<cell-element>',fu)
   enddo

   ! field 06 : a
   write(fu,'(a)')      '# name: a'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(s)
   do is=1,size(s)
     call write_octave(s(is)%a,'<cell-element>',fu)
   enddo

   ! field 07 : ie
   write(fu,'(a)')      '# name: ie'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(s)
   do is=1,size(s)
     call write_octave(s(is)%ie,'r','<cell-element>',fu)
   enddo

 end subroutine write_s1

!-----------------------------------------------------------------------

 subroutine write_e1(e,var_name,fu)
  integer, intent(in) :: fu
  type(t_e1), intent(in) :: e(:)
  character(len=*), intent(in) :: var_name
 
  integer :: ie
  character(len=*), parameter :: &
    this_sub_name = 'write_e'

   write(fu,'(a,a)')    '# name: ',var_name
   write(fu,'(a)')      '# type: struct'
   write(fu,'(a)')      '# length: 7' ! number of fields

   ! field 01 : d
   write(fu,'(a)')      '# name: d'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(e)
   do ie=1,size(e)
     call write_octave(e(ie)%d,'<cell-element>',fu)
   enddo
 
   ! field 02 : i
   write(fu,'(a)')      '# name: i'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(e)
   do ie=1,size(e)
     call write_octave(e(ie)%i,'r','<cell-element>',fu)
   enddo

   ! field 03 : box
   write(fu,'(a)')      '# name: box'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(e)
   do ie=1,size(e)
     call write_octave(e(ie)%box,'<cell-element>',fu)
   enddo

   ! field 04 : bw
   write(fu,'(a)')      '# name: bw'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(e)
   do ie=1,size(e)
     call write_octave(e(ie)%bw,'r','<cell-element>',fu)
   enddo

   ! field 05 : vol
   write(fu,'(a)')      '# name: vol'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(e)
   do ie=1,size(e)
     call write_octave(e(ie)%vol,'<cell-element>',fu)
   enddo

   ! field 06 : si
   write(fu,'(a)')      '# name: si'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(e)
   do ie=1,size(e)
     call write_octave(e(ie)%si,'r','<cell-element>',fu)
   enddo

   ! field 07 : is
   write(fu,'(a)')      '# name: is'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a,i7)')   '# columns: ',size(e)
   do ie=1,size(e)
     call write_octave(e(ie)%is,'<cell-element>',fu)
   enddo

 end subroutine write_e1

!-----------------------------------------------------------------------
 
 !> Read an octave cell array
 !!
 !! The octave array can be either a column or a row; however, it is
 !! assumed to have only one dimension.
 subroutine read_rlist_al(l,var_name,fu,norewind)
  integer, intent(in) :: fu
  type(t_real_lists), allocatable, intent(out) :: l(:)
  character(len=*), intent(in) :: var_name
  logical, intent(in), optional :: norewind
 
  integer :: i, ierr, inrow, incol
  character(len=100) :: message, type
  character(len=*), parameter :: &
    this_sub_name = 'read_rlist_al'
 
   call locate_var(fu,var_name,ierr,norewind)
   if(ierr.ne.0) then
     write(message,'(A,A,A,I3)') &
       'Problems locating "',var_name,'": iostat = ',ierr
     call warning(this_sub_name,this_mod_name,message)
   else
     read(fu, '(A7,A)')   message, type  ! "# type:"," cell"
     read(fu, '(A7,I10)') message, inrow ! "# rows:"
     read(fu,'(A10,I10)') message, incol ! "# columns:"
     allocate( l(inrow*incol) )
     do i=1,size(l)
       call read_octave_al(l(i)%x,"<cell-element>",fu,.true.)
     enddo
   endif
  
 end subroutine read_rlist_al
 
!-----------------------------------------------------------------------

end module mod_tps_phs_grid

